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Abstract 

We numerically investigate margination of white blood cells and demonstrate the 
dependency on a number of conditions including hematocrit, the deformability of the 
cells and the Reynolds number. A detailed mesoscopic hydrodynamic Helfrich-type 
model is derived, validated and used for the simulations to provides a quantitative 
description of the margination of white blood cells. Previous simulation results, 
obtained with less detailed models, could be confirmed, e.g. the largest probability of 
margination of white blood cells at an intermediate range of hematocrit values and a 
decreasing tendency with increasing deformability. The consideration of inertia effects, 
which become of relevance in small vessels, also shows a dependency and leads to less 
pronounced margination of white blood cells with increasing Reynolds number. 


Author Summary 

White blood cells can only properly function if they adhere to the vessel walls. To 
facilitate the adhesion, white blood cells migrate toward the vessel walls in blood flow 
through a process called margination. We model and simulate the process and analyse 
dependencies on a number of conditions including hematocrit, deformability and 
Reynolds number. A detailed mesoscopic hydrodynamic Helfrich-type model is derived, 
validated, and used for the simulations. Besideds the confirmation of previous results, 
obtained with less detailed models, we show a less pronounced margination of white 
blood cells with increasing Reynolds number. Inertia effects, previously neglected, 
become of relevance in small vessels, where the cells may experience Reynolds numbers 
of order unity or higher. These results are also of relevence for the margination of other 
cell types. 


Introduction 

Various experimental and simulation studies of flowing blood have shown that red blood 
cells (RBCs) concentrate in the center of the blood vessel. This can be explained by a 
lift force, arising from cell-wall hydrodynamic interactions, the high deformability of 
RBCs and their nonspherical shapes, see e.g. [^[^. The lift force results in a migration 
of RBCs towards the center of the vessel and a RBC free layer near the wall. 
Differences in size, shape, and deformability are assumed to lead to different lift forces 
and thus a separation of cells with different mechanical properties within the blood 
vessel. White blood cells (WBCs) have a near-spherical shape and are not very 
deformable and thus mechanically different from RBCs. The lift force on WBCs is 
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expected to be much lower than that on RBCs, which suggests, that WBCs may get 
marginated to the RBC free layer near the wall. This effect requires the interaction of 
RBCs and WBCs and is of utmost importance for the functioning of the immune 
system, which requires the adhesion of WBCs to the vessel wall. 

While in principle understood, detailed investigations show a non-trivial dependence 
of WBC margination on various blood flow properties including hematocrit Hf, flow 
rate, vessel geometry, and RBC aggregation e.g. a pronounced margination within 

an intermediate range of Ht ~ 0.2-0.3, and reduced WBC margination for lower and 
higher Hf. Only recently, such behavior could be explained through simulation studies 
in 2D [^. It is argued that for low WBC margination turns out to be weak due to a 
low concentration of RBCs and thus less interaction, while at high Ht WBC 
margination is attenuated due to interactions of marginated WBCs with RBCs near a 
wall, which significantly limit the time WBCs spend near a wall. This argumentation is 
confirmed by 3D simulations in an idealized blood vessel [^[^. We here do not add any 
new physical explanation, but confirm the simulation results by using a different, more 
realistic modeling approach, which is based on a Helfrich-type curvature-elastic 
model with various constraints concerning membrane inextensibility and area 
conservation, cell-cell interaction and incompressible fluid flow of the blood plasma. 
This approach contains all relevant physical phenomena and only measurable 
parameters and thus allows for quantitative predictions. 

The paper is organized as follows. In Methods and Models we first review existing 
modeling approaches, discussing pros and cons, before we introduce in detail the 
Helfrich-type curvature-elastic model. The effect of various parameters on the 
margination of WBCs is discussed in Results and Discussion. The SI further provides 
information on thermodynamic consistency of the model, a detailed numerical approach 
by adaptive finite elements, a parallelization concept and model validation and 
convergence studies. 


Methods and models 

Previous Models 

Previous simulation studies, which have been performed to describe WBC margination, 
are based on strong assumptions. The first simulation approach assumes an 
incompressible Stokes flow, the cells are modeled with a linear elastic membrane and a 
global area constraint is enforced. A boundary integral formulation is used for numerical 
discretization. More recently, a particle-based Lagrangian approach was used [^|^. Here, 
RBCs and WBCs are described by a network model, where the cells are represented 
through triangulated surfaces. Penalty terms are used to ensure global volume and 
global area conservation as well as local area conservation for each surface element. The 
approach thus guarantees inextensibility for sufficiently small surface elements. Each 
membrane point is connected to the fluid through viscous friction. The dynamics of the 
fluid flow is described by the smoothed dissipative particle dynamics (SDPD) method, 
an approximation for the Navier-Stokes equations which is only precise, if the particle 
density is large enough. Furthermore, the incompressibility of the fluid is not 
guaranteed a priori and has to be controlled. In a finite element approach is used for 
the RBCs, which are modeled as biconcave capsules and a Lattice-Boltzmann method 
for the fluid flow. The problems are coupled through an immersed boundary method. 

Our approach will improve on these methods by considering RBCs and WBCs using 
a Helfrich-type curvature elastic model with an inextensibility constraint, which goes 
beyond the linear elastic membrane, network model or biconcave capsule approach 
considered previously. Furthermore we will consider the full Navier-Stokes equations to 
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account also for inertia effects, as RBCs in small vessels may experience Reynolds 
numbers of order unity or higher, especially, if the vessels are constricted due to diseases 
such as thrombosis, e.g. [9 10 . 


Helfrich-type models 

Helfrich-type modeling approaches have been applied to understand the complex 
motions and shape changes RBCs undergo within a flow field, e.g. tank-treating (TT) 
and tumbling (TB) motion [^. Within a low Reynolds number regime, the Stokes limit 
is valid and various numerical approaches have also been considered in this limit to 
analyze the TT and TB motion 12 - 20 . All models consider a membrane free energy 


£ = 


- Ho)^ dr 


( 1 ) 


with membrane r(t), mean curvature spontaneous curvature i^o, and normal 
bending rigidity 6^. Lagrange multipliers are used to enforce a global area constraint or 
the stronger inextensibility constraint. The jump condition for the fluid stress tensor 
S = So,i = —pi + with pressure p, fluid viscosity z/q, cell viscosity Vi and 

deformation tensor D = Vv -|- (Vv)^, with velocity v, along the membrane than reads 


5£ 

[S ■ n]r=^ + XgiobaiHn 

ss 

[S • ^]r=^ + MocaiHn + Vr^ocai 


global area constraint, (2) 

local inextensibility constraint, (3) 


with outer normal n and the surface gradient Vr = PV, where P = I — n (g) n denotes 
the projection operator. The Lagrange multipliers Xgiobai and Xiocai are functionals of 
the fluid velocity v and are obtained by requiring ^ /p dT = Jp iLv • n dT = 0 for the 
global area constraint and Vr * v = 0 along T{t) for the local inextensibility constraint. 
The jump condition for the velocity reads [v]p = 0. 

The linearity of the Stokes problem allows for efficient decoupled algorithms to solve 
for the Lagrange multipliers 18 20 . However, considering blood flow in small blood 
vessels, the Reynolds number at the scale of the RBC, can be of order unity and the 
Stokes limit is at least questionable. Modeling approaches, which consider also inertia 
effects, have recently been introduced 


21- 23 . All have found that the classical TB 


behavior is no longer observed at moderate Reynolds numbers. As such, a suppression 
of TB motion could have far reaching consequences also for the interaction of RBCs and 
thus also the margination of WBCs. We will in the following consider the full 
Navier-Stokes equations, which read inside and outside the cells 


p(9tv + V • Vv) — V • S = 0 
V-v = 0 


( 4 ) 

( 5 ) 


with density p = po 
used by 
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The global area constraint can be treated explicitly, which was 
within a front tracking method, proposed by 25 and used by 


e.g 

for a phase-field model and also considered in for a level-set approach. The local 


inextensibility constraint is more delicate and leads to additional non-linearities, which 
are so far only considered within a level set approach in 21 28 and within a phase field 


approximation in 23 . All these models consider only one cell and it will be the main 
modeling contribution to extent the described approaches in 
interactions of cells in an efficient way. 


23 


to model also 
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Hydrodynamic phase field models 

The method introduces auxiliary phase fields (pi that distinguish the inside and the 
outside of each cell i = 1,..., n. The inside and the outside are separated from each 
other by a diffuse layer, which marks the membrane. The phase field variables are 
defined as 

(piit, x) := tanh ( 6 ) 

where e characterizes the thickness of the diffuse interface and ri{t,x.) denotes the 
signed-distance function between x G O and its nearest point on r^(t) the membrane of 
cell i. Depending on Vi we label the inside with <pi ^ 1 and the outside with pi ^ —1. 
r^(t) is then implicitly defined by the zero level set of pi. The cell phase can thus be 
defined as pceii ~ 1 with pceii = maXxGo(0i, • • •, Pn) fhe fluid phase as 

00 = —Pceii ~ 1 - 

The dynamics are now governed by equations that couple these phase fields to the 
actual physical degrees of freedom. We consider the diffuse nondimensional Helfrich 
energies 


£i{Pi) 


2ReBe, 


/ ■ 

Jn e 


- eA(/)i 






( 7 ) 


with Reynolds number Re = and bending capillary numbers Be^ = ? 

where U denotes a characteristic velocity, L a characteristic length and the bending 
rigidity of cell i. In formal convergence for e ^ 0 to the nondimensional form of the 
sharp interface energy in eq. 0 is shown. 


Instead of a direct extension of the models in 23 which considers an L^-gradient 


flow for pi and enforces a volume and local or global area constraint by Lagrange 
multipliers, we here introduce a gradient flow, which directly ensures volume 
conservation and treat the global area constraint using a penalty approach. Such an 


approach was already considered for one cell, but without flow interactions, in 30 . It 
has the advantage to keep the equations local. We introduce the nondimensional 
penalty energies 




2ReBe. 




( 8 ) 


with penalty parameter c and initial and desired area of cell i, ^o,i and 
A{pi) = ||V0ip + ~ respectively. The last term converges to 
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/p. if e ^ 0 , see 

In addition, we requir^n interaction energy ^i^t, to be defined below, such that the 
overall energy can be written as 

n 

. . . ,(pn) = + £i,area{4>i)) + £int{<t>l, • • • , 4>n)) (9) 

i=l 


and the evolution equations for pi read 

+ V • V(l>i = 7 A(/)§ 


( 10 ) 


with a small positive mobility coefficient 7 and the nondimensional chemical potentials 



d^ii^Pi) dSi^area I ^£int{4^1: • • • 5 0n) 
s<pi ^ 


( 11 ) 
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We obtain 


= p p i’i, i’i = - 4(39^1 + ‘^^0<P - IJ-i = eA^i - -{cff - !)((/)* + Hq). 


5(j). 


ReBe^ 


The penalty terms read 

^^i.area ^ 


pp Ki{Ao,i - A{(pi)), Ki=eA4>i -(</>•-l)</>i. 

KeBe^ e 


( 12 ) 


(13) 


Now, we have to consider the interaction terms. Interaction in principle is 
computationally costly, as it turns the problem into a nonlocal one and requires the 
coupling of all phase field variables ,..., and computations of the distance 
between cells. We here consider only steric interactions and model the short range 
repulsion by a Gaussian potential, which in the sharp interface description reads 


n « 

^i,int(ri,... ,r„) = Wa / WjdT, with Wj(x) = exp 


r 2 (x) 


(14) 


J = 1 


where Wj is an interaction function and describes the influence of cell j on its 
environment. The interaction parameter a > 0 determines the strength of the repulsive 
interaction between cell i and cell j with respect to the evolution of cell i. Using eq. 
the signed distance function Vj can be computed within the diffuse interface region as 


e 1 + 
Vn =-= In 


Vx : \(pj{x)\ < 1 


V 2 1 - 

We thus can approximate the short range interaction function Wj as 

fexp , if |^!.j(x)| < 1 

0 otherwise 


Wn 


(15) 


(16) 


and consider the interaction potential within the phase-field description, which reads in 
nondimensional form 


1 f ^ 

• • • : 0 n) = ^ j / ^ ^ 

3^'i 


WjdQ 


(17) 


with B{(j)i) = — 1)^ being nonzero only within the diffuse interface around T^ and 

In = the interaction number. The interaction energy thus reads 


’)•••’) 0n) — ^ ^ (015 • • • 5 0n) 




and we obtain 


• • • 5 0n) 1 1 ' U ^ A \\ 

- 5^^ -= (-/>») 2 . W'i + 


Rein 


1=1 


1 = 1 


As our approximation of Wi is not differentiable, we define 

10 otherwise. 


(18) 


(19) 


( 20 ) 
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Fig.[T] gives a schematic illustration of the interaction terms. The algorithm 
considers only these cells, for which the diffuse interfaces overlap. All other cells do not 
contribute to the interaction. In addition, the most expensive part, computing the 
distance between cells, has been avoided, as this information is already contained in the 
phase field description of the cells. The approach thus scales with n, the number of cells. 
Similar ideas to model interactions within phase field approaches have been considered 
However, only for the interaction of one cell with a fixed substrate. 


m 
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Figure 1. The red and blue colored parts of are in contact with the interfaces of cell 
j and cell k (dashed contours around those cells), where the signed functions rj and 
can be calculated and thus also the interaction functions Wj and Wk- They do not 
vanish in the overlapping regions. 


The nondimensional Navier-Stokes equation reads 

1 "" 

/9((9tV + V ■ Vv) + Vp - — V ■ (ivD) = 


i=l 

V-v = 0 


( 21 ) 
( 22 ) 

with p = 1 and z/ = 2 °°“ + X)r=i ^ ’^' 2 ^ ■ Different densities could be handled in a 

similar way but are omitted here for simplicity. 

To enforce the local inextensibility constraint, we follow the approach of model B 
in [^. The nonlinear evolution equations for (j)i remain, only the Navier-Stokes 
equation has to be extended and now reads 

p(a*v + v-Vv) + Vp-Tv.(zyD) = ^0fV0i + V-(^^^PA,oeai) (23) 

i=l 

V-v = 0 (24) 

with a Lagrange multiplier Aiocai for which we introduce the additional equation 

(25) 


• (z^LVAzocaz) + : vv = 0 
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Symbol 

Description 

Value 

L 

radius of a perimeter-equivalent circular cell 

5 • lO-^' m 

u 

characteristic velocity 

2.25 • 10-5 m/s 

p 

fluid density 

10 ^ kg/m^ 


dynamic viscosity of the fluid 

10-5 Pas 

l^RBC 

dynamic viscosity of the RBC 

10-5 Pas 

^WBC 

dynamic viscosity of the WBC 

5-10-2 Pas 

bN.RBC 

bending rigidity of the RBC 

2-10-1® J 

bN,WBC 

bending rigidity of the hard WBC 

2-10-1® J 

bN,WBC 

bending rigidity of the soft WBC 

2-10-1® J 

€ 

diffuse interface thickness 

0.04 

7 

regularization parameter 

lo-i" 

a 

repulsion parameter 

8.44 - 10-^ N/m 


Table 1. Mechanical and numerical parameters used in the simulations. Mechanical 
parameters correspond to the considered values in 


with ^ > 0 a parameter independent of e. For e 
r = and P : Vv = Vr • v = 0 near F, which was shown in 


0 we obtain AXiocai = 0 away from 
for n = 1 . 
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In SI Analysis we show thermodynamic consistency of the derived models eq. 
, n and eq. and or eq. ^3 ^ and ^ 


10 and 


lj_ for 2 = 1 ,..., n and eq. and or eq. M and ^ A detailed description of the 
numerical algorithm, which uses a finite element discretization in space and an operator 
splitting approach with a semi-implicit discretization in time is given in SI Numerics. 
The considered parallelization approach can be found in SI Implementation and 
numerical tests are provided in SI Benchmark Computations. 


Results and Discussion 


We study WBC margination for different WBC stiffnesses, different hematocrit values 
and different Reynolds numbers. We consider a blood vessel of thickness 20/im and 
length 40/im with periodic conditions on the in- and outflow. The relatively small 
length results from compromising computational efficiency and physical accuracy and 
has been obtained through detailed investigations on the influence of the periodicity on 
WBC margination. We consider RBCs with perimeter 22/im, area 19.5/im^, bending 
rigidity b^^RBC = 2 • 10“^^ J, viscosity vrbc = 1 * 10“^ Pas. WBCs are initially set to 
be circular with radius 5/im. They have a viscosity i^wBC = 50 • 10“^ Pa s. In order to 
study the influence of the stiffness of the WBCs, we consider three types: soft WBCs 
with bjsf^wBC = 2 • 10“^^ J, hard WBCs with b^^wBC = 2 • 10“^^ J and rigid WBCs. 
The last is considered using a fluid particle dynamics (FPD) approach for the WBCs, 
Please note that the FPD approach does not solve the phase field eq. [Tol it 


see 
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rather updates the cell’s midpoint by the averaged velocity inside the cell. The 
interaction strength is constant between all cell types and reads a = 4.24 • 10“^ N/m. 
For the fluid phase, we consider the viscosity uq = 1 ' 10“^ Pas. We consider a constant 
flow rate, which is realized by applying a time-dependent force term F = 
where Fr denotes the Fronde number. If the current flow rate Q{t) is lower or greater 
than the desired flow rate Qq, we increase the force term by multiplying it with the 
ratio of QolQt- The initial force term can be estimated from its Newtonian value: 
1/Fr(t = 0) = 12Qo/{HfRe), where Hi is the channel height. We pick Qq = 15 for all 
simulations, which implies an averaged velocity of 8.44 • 10 “^ m/s. An overview of all 
used parameters is given in Tab. 

In nondimensional units, the computational domain becomes D = [0,8] x [0,4], with 
periodic boundary conditions in the xi direction. The WBC has the radius 1 and is put 
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at (5,2). RBCs are placed randomly such that they do not overlap. The nondimensional 
numbers read Re= 1.125 • 10“^, = 5.3, Bewsc = 0-53 (hard), Be^^^c = 5.3 

(soft), In= 0.1 and Fr(t = 0) = 4 • 10“^. 

We first vary the deformability of the WBC and keep Ht = 0.293 constant. The 
results are presented in Fig. where the lower left diagram shows the X 2 - coordinate of 
the trajectory of the midpoint of the WBC. After an initial phase, the WBC moves 
towards the wall, but only the rigid WBC can attach to the wall, while the soft WBC 
moves away after a certain time. The lower right diagram shows the probability that 
the midpoint of the cell is within the upper part of the channel with height 0.1. The 
results nicely confirm the findings in [^, that WBC margination is high for rigid cells 
and decreases for softer cells. 




- rigid WBC 

- hard WBC 

-soft WBC 


Figure 2. Simulation snapshot at late time for Ht = 0.293 for rigid, hard and soft 
WBC from left to right (top), X 2 coordinate for the trajectory of the midpoint of the 
WBC (bottom right) and the probability that the midpoint of the WBC is inside a 
defined interval (bottom left). X 2 axis is split into 20 intervals of length 0.1. 


The second test concerns the influence of Hf. We vary the number of RBCs, which 
lead to different values of Ht^ ranging from 0.098 to 0.39. Fig. shows the obtained 
results for a rigid WBC and Fig. [^for a hard one. 

For a rigid WBC, margination can be observed for all considered Hf. However, our 
simulations show a lower tendency to move to the wall for the smallest value of 
Ht = 0.098 and the largest tendency for Ht = 0.195 and Ht = 0.293. For the highest 
value Ht = 0.39, the probability slightly decreases. It seems more likely that due to the 
larger number of RBCs interaction between WBC and RBCs are possible also close to 
the wall, which moves the WBC away from the wall, see t = 10, t = 22 and t = 30. This 
results give evidence for a decreasing WBC margination for high Ht, as also observed 
in [^. In the case of a hard WBC, the cell remains in the center and in contrast to Fig. 

no margination occurs for the lowest Ht considered. Increasing Ht leads to WBC 
margination. However, contact with the wall cannot be achieved. We also do not see 
the tendency for decreasing WBC margination for Ht = 0.39. Further increasing 
^N,WBC or Ht is not possible due to numerical reasons. 

The effect of the Reynolds number on WBC margination is shown in Fig. We 
consider Ht = 0.293 and a rigid WBC. Considering a constant flow rate, we obtain 
WBC margination for Re=1.125 • 10“^, Re=0.05 and Re=l. However, the tendency to 
adhere entirely decreases for higher Reynolds numbers. As Reynolds numbers of order 
unity are realistic for small vessels, especially, if the vessels are constricted due to 
diseases, this effect might have severe influences on the functioning of the immune 
system. The simulation results for such Reynolds numbers indicate, that the tendency of 


































Figure 3. Simulation snapshot at late time for a rigid WBC for Ht = 0.098, 

Ht = 0.195 and Ht = 0.39 from left to right {Ht = 0.293 is shown in Fig. (top). X 2 
coordinate for the trajectory of the midpoint of the WBC (bottom right) and the 
probability that the midpoint of the WBC is inside a defined interval (bottom left). X 2 
axis is split into 20 intervals of length 0.1. 




time t probability P{x2) 


- Ht = 0.098 

- Ht = 0.195 

- Ht = 0.293 

. Ht = 0.390 


Figure 4. Simulation snapshot at late time for a hard WBC for Ht = 0.098, 

Ht = 0.195 and Ht = 0.39 from left to right {Ht = 0.293 is shown in Fig. [^, (top). X 2 
coordinate for the trajectory of the midpoint of the WBC (bottom right) and the 
probability that the midpoint of the WBC is inside a defined interval (bottom left). X 2 
axis is split into 20 intervals of length 0.1. 


RBCs to aggregate in the center of the vessel decreases. This increases the concentration 
of RBCs near the wall and thus leads to a stronger interaction with marginated WBCs, 
which, similar to the situation for large Ht, limits the time WBCs spend near the wall. 


Supporting Information 

SI Video 

Simulation video of WBC margination. The video shows WBC margination in a 
vessel for Ht = 0.293 considering a rigid WBC and adopting the simulation parameters 
Re= 1.125 • 10-^, In= 0.1 and Q = 15. 


9 














































































Re= 1.125 • 10“^ 
Re= 0.05 


Figure 5. Simulation snapshot at late time for a rigid WBC and Hf = 0.293 for 
Re=1.125 • 10“^, Re=0.05 and Re=l from left to right, (top). X 2 coordinate for the 
trajectory of the midpoint of the WBC (bottom right) and the probability that the 
midpoint of the WBC is inside a defined interval (bottom left). X 2 axis is split into 20 
intervals of length 0.1. 


SI Analysis 

Thermodynamic consistency. The proposed system of equations fulfill 
thermodynamic consistency. To show this, we consider the kinetic energy Skin = / 
with constant density p = I and the cell energy E and show that the time derivative is 
less or equal to zero: 


015 • • • 5 0n) — ^kin S — 


/ 


VVt 




(pjdtcpidx 


(26) 


with 


= -V • W(f)i + (27) 

= -(v ■ V)v - Vp+ Tv ■ (z/D) + + V ■ (T^PAzoeaO (28) 

i=l 
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which yields 


^tot(v, 4>i,...,(j)n) = J V ■ ^-(v ■ V)v - Vp + ■ {vT>) + ^ </>! 


V</>i 


+v • ^ + 7A<^5) dx 


partial integration, use V • v = 0 


/ -|^|Vv|2 - 7 y] |V,/.S|2 - (Vv : dx 


use eq. 


25 


i=l 

and partial integration 


= / - ^ I Vv|2 - 7 E I' - dx 


< 0 , 


where we have used the identity v x (V X v) = V(|vp) — (v • V)v from which follows 
that f V • (—V • Vv) = 0. 


S2 Numerics 

Time discretization. In order to discretize in time, we explore an operator splitting 
approach. In an iterative process, we first solve the flow problem and substitute its 
solution into the phase-field equations, which are then solved separately with a parallel 
splitting method. We split the time interval / = [0,T] into equidistant time instants 
0 = to < ti < ... and define the time steps r := t^+i — t^. Of course, adaptive time 
steps may also be used. We define the discrete time derivative := 

where the upper index denotes the time step number and e.g := v(tn) is the value of 
V at time tn- For each system, a semi-implicit time discretization is used, which 
together with an appropriate linearization of the involved non-linear terms leads to a set 
of linear system in each time step. 

Space discretization. We apply the finite element method to discretize in space, 
where a P^/P^ Taylor-Hood element is used for the flow problem, all other quantities 
are discretized in space using P^ elements. In each time step we solve: 

1. the flow problem for and 

+ (v^ . V)v^+^ = (29) 

~ + E + V • (30) 

V • v”+i = 0 (31) 

■ {{^eiifVXlYai) + : Vv”+1 = 0 (32) 

where = v{(j)^) and = I — 
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2. the phase field equations for 0' 


n+l 


i = 1,..., n 


+ v"+i • VC+1 = 7A(/>r^\ 
^r +1 _ 


n+l 


ReBe, 


ReBe, 


(+,i-++))+ 




Rein 






j=i 


J = 1 


/ 


+"+i = A++1 - ^(3(++i)2 + 2no</>r - 1)+ 

M+' = eA++i - - !)(++' + ^o) 


,n+l 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


where we either choose ~ or = V • (We 

again linearize the non-linear terms by a Taylor expansion of order one, e.g. 
((^n+l)2 _ i)^n+l ^ ((^n)2 _ i)^n + (3(^n)2 _ _ ^n). 


S3 Implementation 

Implementation. The fully discretized system of partial differential equations is 
implemented using the adaptive finite element toolbox AMDiS 


35 36 . We use an 


adaptively refined triangular mesh Th with a high resolution along the cell membranes 
to guarantee at least five grid points across the diffuse interface. The criteria to refine 
or coarsen the mesh is purely geometric and related to the phase field variables 0^. Due 
to the use of adaptivity, we need to interpolate the old solution defined on Tf^ onto the 
new mesh • To do this without violating the conservation of cell volume, we solve 
{(j/l'Oid^0^ in every adaption step, with test function 0 and defined 

on and on • We use a multi-mesh strategy 37 to virtually integrate the 


first term on the finest common mesh U ^ which guarantees a constant cell 
volume as long as time steps are appropriately chosen. We require the cell membranes 
not to propagate over a whole element within one time step. With this restriction, all 
numerical tests show that (pi dx is conserved. We conduct a shared memory 
OPENMP parallelization, to solve the phase field evolutions via a parallel splitting 
method. Furthermore, each linear system of equations is solved using the direct 
unsymmetric multifrontal method UMFPACK. 


S4 Benchmark computations 


Benchmark - collision of two cells In 23,28 it has been shown that two 


inextensible cells in an extensional flow do not coalesce, even without an interaction 
potential. If this would also be true in more general cases, it would allow to use one 
phase field variable for all cells and drastically simplify the considered model. We 
therefore consider the interaction of two cells in more detail and demonstrate that the 
considered model is indispensible. 

We set up two elliptical cells at (1.45,0.88) and (3.8,0.875), with axis length 0.5v/2 
in X 2 -direction and 0.5 in xi-direction, where the first cell is placed a little bit higher in 
order to have comparable situations at coalescence, with the left cell always on top. The 
computational domain is [0,5.25] x [0,1.75]. We apply Dirichlet conditions on each 
boundary. To provide a collision, we adopt a space dependent volume force 
F = (O.5(0cezz + l)^/i 5 0)^ added to the right hand side of eq. 21, with fi = 1, if 
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xi < 2.625 and /i = — 1 if xi > 2.625, and choose the Fronde number Fr=10“^. We set 
Re= 0.01, Be= 5, mobility 7 = 10“^ and a viscosity ratio z^cezz/^o = 10- The evolution 
of the cells at time 0.2, 0.4, 0.6 and 2.5 is shown in Fig. Each row considers a 
different modeling approach. We consider six cases: one phase field without (a) and 
with (b) inextensibility constraint, two phase fields without (c) and with (d) 
inextensibility constraint and two phase field and an interaction potential without (e) 
and with (f) inextensibility constraint. The last situation describes the used model in 
the paper. The black lines show the zero level set and the outer interface, e.g. the 
[— 0 . 8 , 0 ] level sets, are colored blue in order to visualize the effects at coalescence. 

The simulations show a strong influence of the inextensibility constraint on the 
dynamics, the cells move more slowly. However, the strongest effect can be seen on 
coalescence. For one phase field variable the cells come into contact, but do not merge, 
as merging would cost energy, due to a violation of the tanh profile. With the 
inextensibility constraint, which should inhibit the cells from touching according 
to [^[^ , the cells come into contact only at two points. However, also in this situation 
the phase field is irreversibly connected. The results of 23 28 thus cannot be 
generalized. However, in our approach the inextensibility condition is only 
asymptotically fulfilled, which seems to avoid an entire adhesion but cannot prevent the 
cells from touching. Overall, in our setting a single phase field is not sufficient to 
simulate more than one cell if contact cannot be avoided. With two phase fields but no 
interaction potential the situation is similar. Here, merging of cells by definition is not 
possible and due to the incompressibility of the fluid, overlapping cells should also not 
occur. However, the simulations show that cells touch each other, adhere and overlap 
slightly. This can be reduced by using the inextensibility constraint but not be 
prevented. Thus, this approach is not sufficient as well. Only with an interaction 
potential, contact or adhesion of cells was not observed, neigther without nor with the 
inextensibility constraint. Only the remaining distance between the cells differs and is 
bigger if the inextensibility constraint is considered. 


Benchmark - RBCs in a bifurcate vessel 

In the following, we consider a more realistic setting, where we put eight RBCs in a 
flow inside a bifurcate vessel. We consider similar parameters as before, 

Re= 1.125 • 10“^, BeRBC = 2.5, Id= 0.01, 7 = 10“^ and a viscosity ratio lyceiilj^o = 10. 
Flow is considered through the force term F = (^,0)^, with Fr= 2.4 • 10“^ leading to a 
maximal velocity of magnitude 10. Fig. shows the results for the same eight cases as 
considered above. Also for this situation only the cases with one phase field for each cell 
and an interaction potential lead to acceptable results. Differences in the dynamics still 
can be observed for the case with and without the inextensibility constraint. 
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